Why use Mizer?
Mizer is a tool that can be used to simulate a dynamic size spectrum in an marine ecosystem, subject to changes through time, such as fishing pressure. Multiple interacting species and fishing gears can be defined allowing for a range of fisheries management strategy scenarios, and their ecosystem impacts, to be tested.
Below we explore some examples of what can be done using a previously calibrated Mizer model.
Species Collapse and Recovery in a Heavily Fished Multispecies System
In this example we will read in a calibrated model for the North Sea, which consists of 12 interacting species and a background resource. All species are predators and prey since they feed according to size-based rules and they encounter each other. Fishing is applied as a logistic curve by length, matched to the maturity length \(l_{50}\) of each species.
More information on how this model was calibrated with historical data can be found in our “mizerHowTo” package (tutorial("HTM3")) and in this blogpost.
First let’s take a look at what the model has used in the historical period (from ICES stock assessments).
The species’ maximum fishing mortality rates (which we call for convenience effort) for Cod and Saithe have declined but are not as low as what would be in line with single-species “Fmsy” (a reference level of fishing deemed sustainable from single-species stock assessements) and the fishing mortality rate for Sprat and Dab is very high.
Then you can take a look at how the modelled species biomasses change through time by running the following:
plotlyBiomass(sim, start_time = 1950, end_time = 2018)# size spectra averaged for the final 5 years
plotlySpectra(sim,time_range = c(2015:2018),power=2)Examining changes in the community relative to an unfished state
To be able to examine any changes in the community relative to an unfished state we can re-set effort = 0 and calculate the steady state.
plotlyBiomass(sim0)We can compare the current size spectra (averaged over 2015-2019) to the unfished size spectra to assess whether there is any evidence of a size-structured trophic cascade due to fishing.
plot_relative_biomass(sim0, sim)Here we can see the effect of the reduction in large sized individuals of heavily fished species on the other sizes and species in the model, relative to the unfished steady state. The abundance of some (but not all) of the smaller to medium sizes of prey are a lot higher when their larger predators are removed (note the logarithmic scale). Sprat looks to be much lower.
We can do the same with Biomass to see when, if any, of the species collapse. For simplicity, we use < 0.1 of B/B_unfished as a proxy for a reference point for population collapse. We can add other reference points to this kind of plot, for example a simple rule of thumb for B_msy could be 0.5*B_unfished.
Projecting future recovery scenarios
We may now wish to explore the potential recovery of the larger species and sizes in the system. To do this we set up another scenario, where the model starts with the last time step of the fished scenario.
The effort values in 2018 were
| Sprat | Sandeel | N.pout | Dab | Herring | Gurnard | Sole | Whiting | Plaice | Haddock | Saithe | Cod |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1.24 | 0.0845412 | 0.311 | 1.278 | 0.191 | 0.3773742 | 0.335 | 0.198 | 0.192 | 0.247 | 0.419 | 0.3200014 |
Let’s start a new simulation that begins with the effort from 2018 and projects forward for 50 years. We will apply a linear reduction in effort for a selected species to a target value (here assumed for simplicity to be 0.2 with effort expressed in terms of the species’ fishing mortality rate for fully selected sizes).
To do this we need to work with the effort array (time x gear) to enable changes in effort through time, for this scenario. Here we have a ‘gear’ for each species, since effort in this case were annual single-species fishing mortality estimates.
Time | Effort |
2019 | 1.2400000 |
2020 | 1.1244444 |
2021 | 1.0088889 |
2022 | 0.8933333 |
2023 | 0.7777778 |
2024 | 0.6622222 |
2025 | 0.5466667 |
2026 | 0.4311111 |
2027 | 0.3155556 |
2028 | 0.2000000 |
2029 | 0.2000000 |
2030 | 0.2000000 |
2031 | 0.2000000 |
2067 | 0.2000000 |
2068 | 0.2000000 |
2069 | 0.2000000 |
We can see that when we reduce fishing on Sprat it increases the biomass of this species but also affects the biomass of other species in the community.
Are any species collapsed still?
Exploration options - break into smaller groups:
Cut and paste the above code chunks and edit it to explore your own fishing scenarios.
Alter some of of the model parameters to see how this affects the scenario outcomes.
Using the unfished sim object examine the effects of adding a species to the model
Further Exploration
Example 2: Set up a fishing scenario
We could examine how this fish community would have responded if there had been a completely different historical development of fishing fleets over time.
Here we will examine how some commonly used community indicators of the ecosystem effects of fishing change in response to this scenario.
We can alter the fishing parameters using a function called gear_params() and by changing the effort input.
Let’s take a look at the fishing parameters. Note that the gears in the above model were already setup to be species-specific.
| species | gear | sel_func | l25 | l50 | catchability | initial_effort |
|---|---|---|---|---|---|---|
| Sprat | Sprat | sigmoid_length | 7.6 | 8.1 | 1 | 1.2953333 |
| Sandeel | Sandeel | sigmoid_length | 9.8 | 11.8 | 1 | 0.0651055 |
| N.pout | N.pout | sigmoid_length | 8.7 | 12.2 | 1 | 0.3138000 |
| Dab | Dab | sigmoid_length | 11.5 | 17.0 | 1 | 0.9780000 |
| Herring | Herring | sigmoid_length | 10.1 | 20.8 | 1 | 0.1815000 |
| Gurnard | Gurnard | sigmoid_length | 19.8 | 29.0 | 1 | 0.4625057 |
| Sole | Sole | sigmoid_length | 16.4 | 25.8 | 1 | 0.3738333 |
| Whiting | Whiting | sigmoid_length | 19.8 | 29.0 | 1 | 0.2426667 |
| Plaice | Plaice | sigmoid_length | 11.5 | 17.0 | 1 | 0.1848333 |
| Haddock | Haddock | sigmoid_length | 19.1 | 24.3 | 1 | 0.3015000 |
| Saithe | Saithe | sigmoid_length | 35.3 | 43.6 | 1 | 0.3930000 |
| Cod | Cod | sigmoid_length | 13.2 | 22.9 | 1 | 0.2666675 |
We can group these species together according to the gears they are caught by. Let’s imagine a big super trawler.
# allocate species to gear types
gear_params(params)$gear <- c("super_trawler")Note that catchability is set to 1. This is because the fishing “effort” was here assumed to be the fishing mortality rate of fully selected sizes (see here setFishing).
The previous effort won’t work with these new gears, as it is gear x time. We only have a single gear now, so this is easier to set up.
params <- setFishing(params, initial_effort = 1)Now let’s run two simulations, one with light fishing effort (effort = 0.5) and one heavely fished (effort = 1.5).
sim_light <- project(params, effort = .5, t_max = 100)
sim_heavy <- project(params, effort = 1.5, t_max = 100)plot_relative_biomass(sim_light,sim_heavy)Large individuals are the most affected by the high fishing effort. This create a trophic cascade where smaller individuals, even if they are affected by the same fishing effort, will increase in abundance due to the release in predation from larger individuals.
Example 3: add a species to the North Sea model
Let’s start with the North Sea model as in Example 1.
comment(params@catchability) = NULL # remove comment, needed to add species
params@species_params$biomass_observed = getBiomass(sim)[nrow(sim@n), ] # set biomass_observed
plotCalibration(sim)plotGrowthCurves(sim, species_panel = TRUE)That is set up and behaving OK. Now let’s add another species and see what happens. We call it ‘Colossal squid’, arbitrarily pick the species parameters, and an observed biomass of 110g.
species_params <- data.frame(
species = "Colossal squid",
w_inf = 3e4,
w_mat = 2e3,
beta = 100,
sigma = 2,
k_vb = 0.2,
a = 0.01,
b = 3
)
new_params = addSpecies(params, species_params)## No h provided for some species, so using f0 and k_vb to calculate it.
new_params@species_params$biomass_observed[13] = 1e10 # set observed biomass
new_params = steady(new_params)## Convergence was achieved in 67.5 years.
sim2 = project(new_params, t_max = 5)
plotlyBiomass(sim2, power = 2)plotBiomassObservedVsModel(sim2)Biomass is way off for squid! Use Gustav’s method in this blog post for retuning parameters.
new_params <- new_params |> calibrateBiomass() |> matchBiomasses() |> steady() |>
calibrateBiomass() |> matchBiomasses() |> steady() |>
calibrateBiomass() |> matchBiomasses() |> steady()## Convergence was achieved in 73.5 years.
## Warning in setBevertonHolt(params): For the following species `erepro` has been
## increased to the smallest possible value: erepro[Colossal squid] = 1.97e-07
## Convergence was achieved in 30 years.
## Warning in setBevertonHolt(params): For the following species `erepro` has been
## increased to the smallest possible value: erepro[Colossal squid] = 1.47e-07
## Convergence was achieved in 4.5 years.
plotBiomassObservedVsModel(new_params)Done now! Look at the steady state size spectrum.
plotlySpectra(new_params, power = 2)And look at the growth curves, compared to von Bertalanffy growth curves.
plotGrowthCurves(new_params, species_panel = TRUE)Have a peek at the reproductive parameters, to check.
| species | biomass_observed | R_max | erepro | |
|---|---|---|---|---|
| Sprat | Sprat | 9.703980e+10 | 2.038528e+10 | 0.0097946 |
| Sandeel | Sandeel | 4.636213e+12 | 7.710521e+13 | 0.0000719 |
| N.pout | N.pout | 3.005304e+11 | 2.371267e+10 | 0.0070738 |
| Dab | Dab | 6.764875e+11 | 2.032365e+10 | 0.0105072 |
| Herring | Herring | 8.989985e+10 | 1.985117e+09 | 0.0114390 |
| Gurnard | Gurnard | 7.067464e+10 | 1.833829e+09 | 0.0002563 |
| Sole | Sole | 3.572392e+10 | 3.360454e+08 | 0.0098435 |
| Whiting | Whiting | 9.066888e+10 | 1.132840e+09 | 0.0033342 |
| Plaice | Plaice | 2.397591e+11 | 8.976962e+09 | 0.0017399 |
| Haddock | Haddock | 2.378042e+11 | 1.304191e+09 | 0.0029266 |
| Saithe | Saithe | 1.318700e+11 | 1.264573e+08 | 0.0636085 |
| Cod | Cod | 5.515711e+11 | 7.844713e+08 | 0.4499833 |
| Colossal squid | Colossal squid | 1.000000e+10 | Inf | 0.0000001 |